Add Folders to Path

%%time
import sys, os

# get current directory
path = os.getcwd()

# get parent directory
parent_directory = os.path.sep.join(path.split(os.path.sep)[:-4])

# add utils folder to current working path
sys.path.append(parent_directory+"/Src/utils")

# add integration folder to current working path
sys.path.append(parent_directory+"/Src/integration")

# add FTLE folder to current working path
sys.path.append(parent_directory+"/Demos/AdvectiveBarriers/FTLE")

CPU times: user 234 µs, sys: 129 µs, total: 363 µs
Wall time: 289 µs

# Suppress numba warnings arising from computation falling back into object mode.
from numba.core.errors import NumbaDeprecationWarning
import warnings

warnings.simplefilter('ignore', category=NumbaDeprecationWarning)

Import Data

Import the AVISO-data from the file 'Agulhas_AVISO.mat' stored in the folder 'Data'.

%%time
import scipy.io as sio

#Import velocity data from file in data-folder
mat_file = sio.loadmat('../../../../Data/Aviso/processed_data/Agulhas_AVISO.mat')

U = mat_file['u'][:,:,::7]
V = mat_file['v'][:,:,::7]
x = mat_file['x']
y = mat_file['y']
time_data = mat_file['t'][:,::7]

CPU times: user 18.3 ms, sys: 13 ms, total: 31.3 ms
Wall time: 30.4 ms

Data/Parameters for Dynamical System

import numpy as np

# number of cores to be used for parallel computing
Ncores = 25

# time resolution of data
dt_data = time_data[0, 1]-time_data[0,0]

# periodic boundary conditions
periodic_x = False
periodic_y = False
periodic = [periodic_x, periodic_y]

# unsteady velocity field
bool_unsteady = True

# defined domain
defined_domain = np.isfinite(U[:,:,0]).astype(int)

## compute meshgrid of dataset
X, Y = np.meshgrid(x, y)

## resolution of meshgrid
dx_data = X[0,1]-X[0,0]
dy_data = Y[1,0]-Y[0,0]

delta = [dx_data, dy_data]

Spatio-Temporal Domain of Dynamical System

%%time
# Initial time (in days)
t0 = 0

# Final time (in days)
tN = 25

# time step-size (in days)
dt = .1

time = np.arange(t0, tN+dt, dt)

# length of time interval (in days)
lenT = tN-t0

# longitudinal and latitudinal boundaries (in degrees)
xmin = -3
xmax = 1
ymin = -32
ymax = -24

# spacing of meshgrid (in degrees)
dx = 0.015
dy = 0.015

x_domain = np.arange(xmin, xmax + dx, dx)
y_domain = np.arange(ymin, ymax + dy, dy)

X_domain, Y_domain = np.meshgrid(x_domain, y_domain)

CPU times: user 892 µs, sys: 1.09 ms, total: 1.98 ms
Wall time: 1.01 ms

Velocity Interpolation

In order to evaluate the velocity field at arbitrary locations and times, we must interpolate the discrete velocity data. The interpolation with respect to time is always linear. The interpolation with respect to space can be chosen to be "cubic" or "linear". In order to favour a smooth velocity field, we interpolate the velocity field in space using a cubic interpolant.

%%time
# Import interpolation function for unsteady flow field
from ipynb.fs.defs.Interpolant import interpolant_unsteady

# Interpolate velocity data using cubic spatial interpolation
Interpolant = interpolant_unsteady(X, Y, U, V, time_data, method = "cubic")

CPU times: user 48.9 ms, sys: 17.2 ms, total: 66.1 ms
Wall time: 65.9 ms

Finite-Time Lyapunov Exponent (FTLE)

Next, we compute the FTLE over the meshgrid over the given time-interval. We iterate over all initial conditions and first calculate the gradient of the flow map using an auxiliary grid. 'aux_grid' specifies the ratio between the auxiliary grid and the original meshgrid. This parameter is generally chosen to be between $ [\dfrac{1}{10}, \dfrac{1}{100}] $. Subsequently, we compute the Cauchy Green strain tensor. From the Cauchy-Green strain tensor we can then compute the FTLE. The iteration over the meshgrid is parallelized.

The FTLE is given by: \begin{equation} \mathrm{FTLE}_{t_0}^{t_N} = \dfrac{1}{t_N-t_0}\log(\sqrt{\lambda_{2}(C_{t_0}^{t_N}(\mathbf{x})}), \end{equation} with $ \lambda_{2}(C_{t_0}^{t_N}(\mathbf{x})) $ denoting the maximum eigenvalue associated to the Cauchy Green strain tensor over the time-interval $ [t_0, t_N] $

%%time
# Import gradient of flow map
from ipynb.fs.defs.gradient_flowmap import gradient_flowmap

# Import package which computes eigenvalues/eigenvectors
from ipynb.fs.defs.eigen import eigen

# Import package for progress bar
from tqdm.notebook import tqdm

# Import package for parallel computing
from joblib import Parallel, delayed

# Import package for computing CauchyGreen (CG) strain tensor
from ipynb.fs.defs.CauchyGreen import CauchyGreen

# Import package for computing Finite-Time Lyapunov Exponent (FTLE)
from ipynb.fs.defs.FTLE import _FTLE

# Define ratio of auxiliary grid spacing vs original grid_spacing
aux_grid_ratio = .1 # [1/10, 1/100]
aux_grid = [np.around(aux_grid_ratio*(X_domain[0, 1]-X_domain[0, 0]), 5), np.around(aux_grid_ratio*(Y_domain[1, 0]-Y_domain[0, 0]), 5)]

def parallel_FTLE(i):
    
    FTLE_parallel = X_domain[0,:].copy()*np.nan
    
    for j in range(X_domain.shape[1]):
        
        # set initial condition
        x = np.array([X_domain[i, j], Y_domain[i, j]])
    
        # compute gradient of flow map from finite differencing 
        gradFmap = gradient_flowmap(time, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta, aux_grid)
        
        # gradFmap has shape (2, 2, len(time)) --> we need gradient of flow map from t0 to tN
        gradFmap_t0_tN = gradFmap[:,:,-1]
        
        # compute CG
        C = CauchyGreen(gradFmap_t0_tN)
        
        FTLE_parallel[j] = _FTLE(C, lenT)
        
    return FTLE_parallel

FTLE = np.array(Parallel(n_jobs=Ncores, verbose = 0)(delayed(parallel_FTLE)(i) for i in tqdm(range(X_domain.shape[0]))))

CPU times: user 7.82 s, sys: 1.11 s, total: 8.93 s
Wall time: 10min 58s

######################## PLOT RESULTS ########################
import matplotlib.pyplot as plt

# Figure/Axes
fig = plt.figure(figsize=(16, 10), dpi = 600)
ax = plt.axes()

# Contourplot of TSE over meshgrid of initial conditions
cax = ax.contourf(X_domain, Y_domain, FTLE, cmap = "gist_gray", levels = 600)

# Axis Labels
ax.set_xlabel("long (°)", fontsize = 16)
ax.set_ylabel("lat (°)", fontsize = 16)

# Ticks
ax.set_xticks(np.arange(np.min(X_domain), np.max(X_domain), 1))
ax.set_yticks(np.arange(np.min(Y_domain), np.max(Y_domain), 1))

# Colorbar
cbar = fig.colorbar(cax, ticks = np.linspace(0, .4, 9))
cbar.ax.set_ylabel(r'$(\dfrac{1}{d})$', rotation = 0, labelpad = 10, fontsize = 10)

# Title
ax.set_title(r'$ \mathrm{FTLE}$'+f'$_{{{int(time[0])}}}^{{{int(time[-1])}}}$'+r'$(\mathbf{x})$', fontsize = 20)

plt.show()